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Abstract 



We apply a functional-integral formalism for Markovian birth and death 
processes to determine asymptotic corrections to mean-field theory in the 
Malthus-Verhulst process (MVP). Expanding about the stationary mean- 
field solution, we identify an expansion parameter that is small in the limit 
of large mean population, and derive a diagrammatic expansion in powers 
of this parameter. The series is evaluated to fifth order using computational 
enumeration of diagrams. Although the MVP has no stationary state, we 
obtain good agreement with the associated quasi- stationary values for the 
moments of the population size, provided the mean population size is not 
small. We compare our results with those of van Kampen's O-expansion, 
and apply our method to the MVP with input, for which a stationary state 
does exist. 
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I. INTRODUCTION 



The need to analyze Markov processes described by a master equation arises frequently 
in physics and related fields [1,2]. Since such equations do not in general admit an exact 
solution, approximation methods are of interest. A widely applied approximation scheme 
is van Kampen's '^-expansion', which furnishes corrections to the (deterministic) mean- 
field or macroscopic description in the limit of large effective system size [1]. Another 
approximation method, based on a path-integral representation for birth-and-death type 
processes, was proposed by Peliti [3] . This approach was recently reviewed and extended 
[4], and applied to derive a series expansion for the activity in a stochastic sandpile [5]. 

It should be noted that while effectively exact results can be obtained via numeri- 
cal analysis of the master equation, the calculations become extremely cumbersome for 
large populations or multivariate processes. A further limitation of numerical analyses 
is that they do not furnish algebraic expressions that may be required in theoretical de- 
velopments. For these reasons it is highly desirable to study approximation methods for 
stochastic processes. 

In the present work we apply the path-integral based perturbation approach to a 
simpler problem, namely, the Malthus-Verhulst process (MVP), a birth-and-death process 
in which unlimited population growth is prevented by a saturation effect. (The death rate 
per individual grows linearly with population size.) This is an important, though highly 
simplified model in population dynamics. Although the master equation for this process is 
readily solved numerically, the model serves as a useful testing ground for approximation 
methods. A lattice of coupled MVPs exhibits (in the infinite-size limit) a phase transition 
belonging to the directed percolation universality class. 

In the perturbation approach developed here [3,4], moments of the population size 
n are expressed as functional integrals over a pair of functions, ip(t) and ip(t), involving 
an effective action. The latter, obtained from an exact mapping of the original Markov 
process, generally includes a part that is bilinear in the functions ip{t) and ip(t), whose 
moments can be determined exactly, and 'interaction' terms of higher order, that are 
treated in an approximate manner. In the present approach, the interaction terms are 
analyzed in a perturbative fashion, leading to a diagrammatic series. With increasing 
order, the number of diagrams grows explosively, so that it becomes convenient to devise 
a computational algorithm for their enumeration and evaluation. Elaboration of such 
an algorithm does not, however, require any very sophisticated techniques, and could in 
fact be applied to a variety of problems. This is the approach that was applied to the 
stochastic sandpile in Ref. [5]. In the latter case, the evaluation of diagrams involves 
calculating multidimensional wave-vector integrals. The present example is free of this 
complication, allowing us to derive a slightly longer series than for the sandpile. 

In this work we focus on stationary moments of the MVP. The series expressions are 
apparently divergent, but nevertheless provide nearly perfect predictions away from the 
small-population regime, as compared with direct numerical evaluation of quasi-stationary 
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properties. One might suppose that the divergent nature of the perturbation series is due 
to the MVP not possessing a true stationary state. (The process must eventually become 
trapped in the absorbing state, although the lifetime grows exponentially with the mean 
population [6].) Applying our method to the MVP with a steady input, which does possess 
a stationary state, we find however that the perturbation series continues to be divergent, 
although again providing excellent predictions over most of parameter space. 

In the following section we define the Malthus-Verhlst process, explain the perturbation 
method and report the series coefficients for the first four moments, up to fifth order in 
the expansion parameter. In Sec. Ill we briefly compare these results to those of the 
f2-expansion. Then in Sec. IV we present numerical comparisons of our method (and of 
the fi-expansions) against quasi- stationary properties. We apply our method to the MVP 
with input in Sec. V, and summarize our findings in Sec. VI. 



II. PERTURBATION THEORY FOR THE MALTHUS-VERHULST PROCESS 

Consider the Malthus-Verhulst process (MVP) n(t), in which each individual has a 
rate A to reproduce, and a rate of \i + u(n — 1) to die, if the total population is n. By 
an appropriate choice of time scale we can eliminate one of these parameters; we choose 
to set /j, — 1. Then in what follows we use a dimensionless time variable t' = fit and 
dimensionless rates A' = A//i and v' = vj ji. From here on we drop the primes. 

The mean-field or rate equation description of the process is 

x — (A — l)x — vx 2 (1) 

where x = (n(t)), with the nontrivial stationary solution x = (A — for A > 1. We 
are interested in deriving systematic corrections to this result, and in calculating higher 
moments of the process. 

Our starting point is the expression for the r-th factorial moment of a general process 
taking non-negative integer values, 

(n r (t)) f = e-^ r \C=P), (2) 

where, for simplicity, we have assumed an initial Poisson distribution with parameter p 
(see Eq. (108) of Ref. [4]), so that the probability generating function at time zero is 
$o(-2) — e^ 2-1 ). Here, the kernel is given by the functional integral, 

Ut\0 = ip^^) = j^j Vi>m r mi>l=ieM-Si\ , (3) 



(equivalent to Eq. (106) of [4]), with 
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F[i),i)\ z= i = exp 



- / dt'^[d t/ + (l-\)}^+( 



(4) 



containing the bilinear part of the action, and, in the case of the MVP, the "interaction" 
part, 



S I= I dt'\-\i) 2 i)+ui)(l+i))i) 2 }= I dt'£i{t'), 
'o Jo 



(5) 



(see Eq. (54) of [4]). 

Our goal is to expand each factorial moment about its mean-field value. To this end, 
consider the shift of variable, 



ip(t) =n + (j>(t) 



(6) 



where n is a real constant. On performing this shift, the argument of the exponential 
(i.e., of the factors T\\\)^\ z= \ and exp[— Si]) in Eq. (3) becomes, 



C+ / dt'{-^[d v + l-\ + 2un]<p-n{l-\ + un)^ 
Jo 

+ n(X - vn)i) 2 + (A - 2vn)^ 2 <p - v^<p 2 - v^ 2 (j) 2 } 



(7) 



This simplifies if we let n = (A — 1)/V, which eliminates the term oc ip. Introducing 
w = A — 1 (equal to — w as defined in [4]), the argument of the exponential is: 



C + / dt'{-^[d t ,+w](j) + n^ 2 



+ (2-A)^ 2 - - z/^ 2 2 } 



(8) 



We recognize the exponential of ( plus the first term in the integrand as Fty, 4>] z =i] the 
remaining terms then represent — S' T , the new effective interaction, which will be treated 
perturbatively. The four terms of Sj are represented graphically, following the conventions 
of Ref. [4], in Fig. 1. 
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Fig.l. Vertices in the perturbation series for the 




We refer to these vertices as source, bifurcation; conjunction and 4-veifex, respectively. 



A. Perturbation expansion 




The analysis of the MVP now follows the lines of [4]. Let 

[A} = [nf>[v4,A(<f>J)F[<j>,4>}, 



(9) 



denote the free expectation of any function A of and ip. From the discussion of Sec. 4 
Ref. [4] we have, 



m) r ] = (P - n) r e 



r „—rwt 



(10) 



(here we used 0(0) = p — n, and have already canceled the factor in U® with the 
corresponding factor in the inital generating function, <& (C))> 



(ii) 



and 



[0(*i)A*2)] = e(ti-t 2 ) e - 



-U>(tl-*2) 



(12) 



Consider now the series for the mean population size (n(t)). The zeroth-order term is 
simply 



(n(t)) = n + [<p(t)] —n + (p — n)e 



(13) 



which converges to the mean-field result as t — > oo. Similarly, the leading term in {n r (t)) f 
is rf, so that the stationary probability distribution, in mean-field approximation, is 
Poisson with parameter n. 
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Corrections involving Sj are conveniently represented as diagrams with all lines leaving 
vertices contracted with ingoing lines, either at vertices to the left, or at the "sink" lying 
to the left of all vertices. (In the calculation of (n r (t)) f the sink is a point with r incoming 
lines.) The lowest-order diagram in (n(t)) is depicted in Fig. 2. 




Fig. 2. One- and two-vertex diagrams in the series for (n(t)). 

This diagram makes the following contribution to (n(t)): 

-v f dt ie - w ^\p - n) 2 e~ 2wtl = --(p - n) 2 e~ wt {l - e~ wt ). (14) 
Jo w 

The four diagrams (in the series for (n(t))) having two vertices are also shown in Fig. 
2. Only the first of these four makes a nonzero contribution to (n^) = lim t _ voo (n(i)), 
namely, 

-2un f dh f 1 dhe-^-^e- 2 ™^-^ = --(1 - e~ wt ) 2 . (15) 
Jo Jo w 

(The combinatorial factor 2 represents the number of ways the lines may be contracted 
between source and bifurcation.) 

If we are only interested in stationary properties, it is convenient to use the Laplace 
transform. Let fjj be the n-fold integral over time variables in a given n-vertex diagram 
D. Noting that time-dependent factors are associated with each propagator and each 
uncontracted incoming line, we see that fn is of the form: 

ft rh rt n -i 

f D (t)= dh dt 2 ... dtne -ai(t-t 1 )-a 2 (t 1 -t2)...-c^(t n - 1 -t n )-p 1 t 1 -----hu ( 16 ) 
Jo Jo Jo 

where fa is w times the number of uncontracted lines incident on vertex i. (In the first 
graph of Fig. 2, [3\ = 2w.) The factors are given by w times the number of lines 
(propagators) running between vertices i and i — 1 (with i — representing the sink), 
regardless of where these lines originate or terminate. Now consider 

/•OO 

f D (s) = / dte' st f D (t). (17) 
Jo 
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Using Eq. (16) we may write 

fo(s)= dt dti... dt n exp[-(a 1 + s)(t-ti) - (a 2 + Pi + s)(t 1 -t 2 ) 

Jt! Jt 2 JO 

-(a 3 +p 1 +p 2 + s)(t 2 -t 3 ) (a„+/3i + - ••+/3„-i+s)(* n _i-* n ) 

-(/?! + • • • + /3 n + S)t n ] 

= [(a 1 +s)(a 2 +(3 1 +s) ■ ■ ■ + - • ■+(3 n - 1 +s)((3 1 + - ■ ■+(3 n +s)}- 1 . (18) 

The contribution of D to (n^) is proportional to 

f D = hmf D (t) = ^msf D (s), (19) 

by the Final Value Theorem of Laplace transforms. This is zero unless = 0, Vi = 1, n. 
In the latter case, 

J D = [(«! + s)(a 2 + s) • • • (a n + s)]" 1 . (20) 

Thus the only diagrams contributing to (n^) (and by extension, to (ri^)) are those free 
of uncontracted lines. This is in fact evident on physical grounds: each such line carries a 
factor p — n, whereas stationary properties cannot depend on the initial mean population 
p. In diagrams contributing to (n^), the first vertex (i.e., immediately to the right of the 
sink) must be a conjunction, while the n-the vertex must be a source. 



B. Diagrammatic analysis for the stationary mean population 

The contribution of a diagram D to (rioo) is the product of three factors: f D discussed 
above; a combinatorial factor counting the number of contractions consistent with the 
diagram topology; the product of vertex- associated factors shown in Fig. I. 

Certain infinite classes of diagrams may be summed up exactly. Consider the sequence 
shown in Fig. 3: between the source and conjunction we insert any number of 4- vertices 
or bifurcation-conjunction pairs. 




Fig. 3. A summable set of diagrams in the series for (n(t)). 
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A short calculation shows that the contribution of this series is given by: 



OO I- 



( 



2-A 



) 



1 



w 



-y (--) 

n=0 L 



1 + 



W 



w + vjw 



(21) 



In this way, by multiplying the n = 2 contribution, — — , by k = (1 + u/w 2 ) -1 , we have 
included all diagrams of the form of Fig. 3, i.e., diagrams in which all lines exiting vertex 
% are contracted on vertex i — Vi = 1, ...,n. With this simple replacement such reducible 
diagrams are no longer to be included explicitly. 

A similar observation allows us to add arbitrary sequences of 4-vertices or bifurcation- 
conjunction pairs to any diagram of four or more vertices. (Diagrams outside the series 
depicted in Fig. 3 have n > 4 vertices.) The diagram consists of a 'body' containing 
vertices 2, ...,n — 1 linked to a conjunction (vertex 1) and a source (vertex n). Call the 
factor associated with the body A, so that the diagram makes a contribution of —uAn to 
(rioo). A family of diagrams may be constructed by adding sequences, as above, between 
vertex n and the body; the sum of these contributions is 



By the same reasoning we may insert arbitrary sequences between the first vertex and the 
body, yielding the same set of contributions. As a result, we may multiply the contribution 
of any diagram (not included in the sequence of Fig. 3) by k 2 , and thereby automatically 
include all diagrams with the same body but arbitrary linear sequences between the first 
vertex and the body, and between vertex n and the body. Such diagrams are no longer 
included explicitly. This procedure will be called "dressing" the conjunction (vertex 1) 
and the source (vertex n). 

At this stage we have essentially exhausted the simplifications (via "dressing" or sum- 
ming sets of trivial variations to a diagram) that can be realized in the Laplace transform 
representation. Summarizing, the contributions to {n^ are as follows. 

1. The mean-field contribution n. 

2. The set of diagrams shown in Fig. 3, giving —k/w. 

3. Diagrams of four or more vertices, with no dangling lines, and subject to the restrictions 
mentioned above. Each such contribution is to be multiplied by k 2 . 

Including diagrams of up to four vertices we find: 





(rioo) =n 1 j 



kv 2k 2 v 2 




H 



w 




(23) 
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where the higher order terms come from diagrams of five or more vertices. This result 
suggests that we adopt 



e = 1 - k = -— (24) 

W 2 + V 

as the expansion parameter. Up to diagrams of four vertices we have 

(rO =n [l-e + 2(2- A)e 2 + • • •] . (25) 

The first correction —lit oc 1/A, as A — > oo. 

Note that z/ oc e to lowest order. The lowest order in e at which a given diagram 
contributes to (n^), when expressed in the form of Eq. (25), may be found as follows. 
First observe that the diagram carries a factor v c+ f~ s+1 where c, / and s represent the 
number of conjunctions, four-vertices and sources, respectively, and the contribution of 1 
in the exponent is due to the prefactor n. Certain constraints exist between the numbers 
of vertices. Equating the total number of lines emanating from vertices with the total 
number entering, we find 

2(s + f + b)+c = 2(f + c) + b + l (26) 



where b is the number of bifurcations and the 1 on the r.h.s. represents the sink. Thus 
c = 2s + b — 1. On the other hand the total number of vertices is n = s + b + / + c, and 
using this we find the power of —v tobec + /- s + l = n- c = m. 

A simple analysis yields the algebraic factor associated with a given diagram having 
n > 4 vertices: 

Fai 9 = (-1) C+ V(2 - A) V" n -V- C . (27) 
This is readily expressed in terms of the parameter e as: 

TO 

F al9 = (-1)^(2 - A) V ?r= ^ =5 . (28) 



Since we intend to organize the expansion in powers of e, it is of interest to know 
which values of n correspond to a given m. We begin by noting that in diagrams of four 
or more vertices, the second vertex (from the left) must be a conjunction, just as the 
first is. (If it were a bifurcation or a four-vertex the result would be a diagram already 
included by dressing the first vertex.) Thus c > 2, and there are in fact diagrams with 
just two conjunctions, for any n > 4. We therefore have m<n-2orn>m + 2. To 
find an upper bound on n, we find an upper bound on c. Recall that c = 2s + b — 1. To 
maximize c we let b = f = (a diagram consisting of only sources and conjunctions), in 
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which case c+ s = n, yielding m = n — c = s = (n + l)/3, or n < 3m — 1. Thus for m = 2 
we have contributions from diagrams of 4 and 5 vertices, while for m = 5 the number of 
vertices ranges from 7 to 14. 

In order to evaluate the contributions from diagrams of five or more vertices we have 
developed an enumeration code. Reducible diagrams are eliminated by imposing the 
following rules: 

1) The second vertex, like the first, is a conjunction. 

2) Vertex n — 1 cannot be a 4-vertex. 

3) If vertex n — 1 is a conjunction, vertex n — 2 must be a source. 

Using the enumeration code we are able to evaluate contributions up to C(e 5 ). At this 
order there are approximately 8 x 10 s diagrams. (There are 9, 1317 and 594 339 diagrams 
at orders 2, 3 and 4, respectively.) This explosive growth in the number of terms prevents 
our going to higher order. The fifth-order calculation requires about two days on a fast 
PC. 

Writing Eq. (25) in the form (noo) = n[l + J2 m >i9m£ m ], and using u = 2 — A, we 
have: 

9i = -1, 
g 2 = 2u- 5, 

g 3 = -10m 2 + (55 - Aw)u + 19w - 60, 

g A = 78m 3 - (689 - 62w)u 2 + (1665 - 50Aw + 8w 2 )u 
-65w 2 + 745w- 1165, 

g 5 = -750m 4 + (9437 - 880w)u 3 - (37078 - 10259w + 266w 2 )u 2 
+ (57620 - 32328m; + 3117m; 2 - 16m> 3 )m 

+211m; 3 - 6040m; 2 + 27786m; - 29390 (29) 



A particularly simple case is A = 2, for which we have 

(n) = -[1 - e - 5e 2 - 41e 3 - 485e 4 - 7443e 5 - ...] (30) 

The ratios of successive coefficients in the series are: 1, 5, 8.2, 11.829 and 15.346, sug- 
gesting unlimited growth and therefore a divergent series. Numerical results (see below) 
nevertheless reveal excellent agreement with quasi-stationary properties for A sufficiently 
large. 



10 



C. Higher moments 



We turn now to the expansion of higher stationary factorial moments. From Eqs. (3) 
and (6) we have 



(n r )f = [(n + <l>re- 8 ']. (31) 



Expanding the product, we have first the mean-field contribution ff, and then a series of 
terms involving diagrams. The term oc [<p q e~ Sl } involves (for q < r) diagrams that already 
appeared in the series for the q-th factorial moment. Thus for r = 2 we have 

(n 2 ) f = n 2 + 2n[<f)e- Sl ] + [(p 2 e~ Sl }. (32) 



Note that [0e~ 5/ ] is simply T>i, the sum of all diagrams with a single line entering the sink, 
that is, the diagrammatic series for (n). Consider now the final term in Eq. (32), the series 
V 2 of diagrams with two lines entering the sink. Recalling that, in all diagrams in Vi, the 
first (leftmost) vertex is a conjunction, we see that there is a one-to-one correspondence 
between T>\ and T>2- For a given diagram, making a contribution of A to T>i, there is 
a corresponding diagram (with the conjunction and one-line sink replaced by a two-line 
sink) that contributes —wA/u = —nA to V 2 , that is, V 2 = —nD\. Thus we have 

(n 2 ) f = n 2 + nV 1 

= n(n), (33) 

where in the second line we used (n) = n + Vi. Using Eqs. (29) and (33), we find that 

var(n) - (n) = n 2 e[l - (2u - 4)e + • • •]. (34) 



For a Poisson distibrubtion the difference is zero; here it approaches \ jv as A — > oo. 

For r > 3 we do not have a simple relation between D r and D±, so the diagrammatic 
series must be evaluated to the desired order. In general we have 



3=0 V ^ ' 



(35) 



where T>\ 



~ Sl ], with T>q = 1. (Note that for r = 3 the j = I and j '• = 2 terms 



cancel.) As in the case of (n), we include an overall factor of rf, and write 



n 



E 



U\ 3 

w/ 



(36) 
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Note that for r > 3 we can dress the leftmost source, as before, but that it is no longer 
possible to dress the sink. Thus the algebraic factor associated with an arbitrary diagram 
in (v/w) r V r is 




ku w 



(37) 



An analysis along the lines presented in the r = 1 case leads to to = n — c (as before) for 
the order in u, and to s — n — r — f, so that 



The limits on n, for fixed order to, are to < n < 3to — r. (Note that n = to is possible 
for r > 3 because there are diagrams with no conjunctions.) 

As in the case r = 1, we have constructed an enumeration code to evaluate the 
corrections due to diagrams. Writing the contribution to (n r ) / coming from T> r as 
^ E m >i 9r,m£ m , we have, to r = 4 and to = 5: 



03,1 = 0, 

5(3,2 = 2u - 5, 

53 3 = -10m 2 + (57 - Aw)u + 19w - 65, 

g 3A = 78m 3 - (699 - 62w)u 2 + (1722 - 508u> + 8w 2 )u 
-65w 2 + 764w; - 1230, 

g 3 5 = -750m 4 + (9515 - 880w)m 3 - (37777 - 10581w + 266w 2 )u 2 
+ (58432 - 32836w + 3125m; 2 - 16m; 3 )m 

-30560 + 28550m;- 6105m; 2 + 211m; 3 (39) 



91.1 = 0, 

94.2 = 3, 

54.3 = 6m 2 — 41m — 9m; + 53, 

g A A = -54m 3 + (547 - 30m;)m 2 - (1452 - 330m;)m 
+27m; 2 - 576m; + 1088, 

g A b = 570 M 4 - (7785 - 564m;)m 3 + (32515 - 7799m; + 114m; 2 )m 2 
-(52118 - 26354m; + 1895m; 2 )m 

+28058 - 24234m; + 4361m; 2 - 81m; 3 . (40) 



Falg = (-l) C+f U b W 



(38) 



m— 1 ' 



and 
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Numerical comparisons of these series against quasi- stationary properties of the MVP will 
be discussed in Sec. 4. 



III. COMPARISON WITH THE ^-EXPANSION 



A well known method for obtaining approximate solutions to the master equation is van 
Kampen's fi-expansion [1]. The method depends on the system size (or the expected value 
of the stochastic variable of interest) being large, so that fluctuations are small compared 
to the value n(t) predicted by the macroscopic equation. The stochastic variable n is then 
written in the form 

n = n((t) + n 1 / 2 at) (4i) 



where the first, deterministic, term represents the solution to the macroscopic equation, 
with Q denoting the size of the system. The stochastic contribution, represented by the 
second term, is expected on general grounds to scale as fi 1 / 2 . As shown in Ref. [1], ((t) 
satisfies the macroscopic equation, while the probability density IT(£, t) satisfies, to lowest 
order in fi -1 / 2 , a linear Fokker-Planck equation. In the present case the macroscopic 
equation is 

f t = (A - 1)C - K 2 , (42) 



i.e., the Malthus-Verhulst equation, with nontrivial stationary solution £ = (A — \)jv — n. 
The equation governing the evolution of II is, in the present instance, 

^ = (1 + 2K - A)|(^n) + |C(1 + < + A)|5 (43) 

This implies that in the stationary state, £ is Gaussian with mean zero and variance X/u, 
so that, setting the formal expansion parameter Q to unity in Eq. (41), n is Gaussian 
with mean n and variance Xjv. In the perturbation expansion developed in the preceding 
section, n is, to zeroth order, a Poissonian random variable with mean n. Conceptually, 
a Poisson distribution seems preferable to a Gaussian as the reference distribution (since 
n is discrete and cannot assume negative values), but in practical terms we expect the 
difference to be very small. 

To first order in e, we found (n) = n(l — e) and var(n) = n(l + fie). Noting that 
x A-l /A-l\ 2 v X 1 

^ 1 + ^ = — + {—) ( a-i)» + „ = "(a3iP (44) 



we see that the difference from the f2-expansion result, to first order in f2~ 1//2 , is small 
for A > 1. In the f2-expansion, corrections to the moments (£ r ) can be obtained via the 
procedure detailed in Ref. [1]. In the present case one finds 
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(45) 



so that 



(n) = n 



(46) 



The series prediction is 

(n) = n(l — e) = n 



(A-l) : 



+ 0(e 2 ) 



(47) 



so that the two methods agree to first order in e. Extending the f2-expansion to include 
terms of order Q~ l in £, one finds 



(n) = n 



v 



(A -I) 2 (A-l) 4 



(A 2 - 3A + 4) 



(48 



These results are compared against the e series in the following section. 



IV. NUMERICAL COMPARISONS 



In this section we compare the e-series predictions with exact (numerical) results for 
quasi- stationary (QS) properties of the MVP. The latter are obtained via recursion rela- 
tions leading to the QS distribution as detailed in [6]. QS properties are those obtaining 
at arbitrarily long times, conditioned on survival (i.e., the process has never visited the 
absorbing state). Although a condition on survival is not involved in the perturbative 
analysis of Sec. 2, it seems reasonable to compare its predictions with QS properties, 
since the true stationary properties are the trivial ones of the absorbing state (population 
zero, no fluctuations). We note that, as suggested above, the series in e appears to be 
divergent for any set of parameter values, as reflected in the ratios between coefficients of 
successive terms, whose ratios appear to grow without limit. For suitably large values of 
A the ratios are very small, although increasing with order m, so that the series, truncated 
at fifth order, appears to have "converged". The numerical evidence (limited to m < 5, 
of course) points nevertheless to a divergent series. 

In Fig. 4 we compare the QS mean population size with (rioo) as predicted by the 
series, as a function of A, with v fixed at 0.01. For A less than about 1.2, the series is 
useless (it yields a negative population size!), and is inferior to mean-field theory. For 
A > 1.4 on the other hand, we find good agreement with the QS result. For a more 
detailed comparison we plot, in Fig. 5, the difference An between (n) (as given by the QS 
distribution and the e series) and the mean-field prediction n; there is perfect agreement 
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for A > 1.5. (The correction to the mean population size decays oc 1/A for large A. This 
is readily seen from Eq. (23) if we note that ne ~ 1/A for A ^> 1.) Figures 6 and 7 show 
that similar behavior obtains for v = 0.1 and 0.001. The smaller u, the smaller the value 
of A required for agreement between series and QS values. (For v = 0.1 agreement sets 
in for e < 0.024, approximately; in the other cases e < 0.06 is sufficient.) At the point 
where the series and QS results begin to agree, the mean population is not particularly 
large; for v — 0.01, A = 1.4 corresponds to (n) ~ 35. 




1.0 1.2 1.4 1.6 1.8 2.0 

X 

Fig. 4. Stationary mean population size versus birth rate A in the MVP with v = 0.01. Solid 
line: exact QS value; dotted line: mean-field prediction n; dashed line: series to order 0(e 5 ). 
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Fig. 5. Difference An between the stationary mean population and the mean-field value in 
the MVP with v = 0.01. Solid line: exact QS value; dashed line: fifth-order series. 
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Fig. 6. Same as Fig. 5 but for v = 0.1. 



In order to gauge the importance of successive terms, we plot (Fig. 8) An (for v = 0.01) 
as predicted by the series to order e m , for m = 1,..., 5. For A > 1.5 very good agreement 
with the QS mean population is obtained using the result to 0(e 3 ). For this value of A 
the absolute difference between the three- and five-term series is about 0.04 The relative 
difference is about one part in a thousand; the difference rapidly decreases for larger A. 
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Fig. 7. Same as Fig. 5 but for v = 0.001. 
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Fig. 8. Stationary mean population size versus birth rate A in the MVP with v = 0.01. The 
curve exhibiting a minimum near A = 1.3 represents the exact QS value. The other curves 
represent (in decreasing magnitude) the series truncated at first, second,..., fifth order. 
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Fig. 9. Difference Avar(n) between the stationary variance and its mean-field value in the 
MVP with v = 0.01. Solid line: exact QS value; dashed line: fifth-order series. 
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In light of the poor performance of the e series in the vicinity of A = 1, it is natural to 
apply a resummation technique such as Pade approximants. We therefore constructed the 
[2,3], [3,2] and [4,1] approximants to the series 1 + g±e + • — hffee 5 and to the series for the 
logarithm of this expression. None of the Pade approximants yielded any improvement 
over the original series; the approximants are ill-behaved (large and negative) near A = 1 
(they typically exhibit a pole in this region) and reproduce the excellent agreement with 
the QS results whenever the original series does. 
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Fig. 10. As in Fig. 9, but for the third central moment. 

The trends observed for the mean population continue when we compare series and 
QS predictions for higher moments. In Fig. 9 we plot the difference, Avar(n), between 
the variance furnished by these predictions and that given by mean-field theory. (Since 
the latter yields a Poisson distribution, the variance is simply equal to n in mean-field 
approximation.) For A ~ 1 the series prediction is useless; close agreement with the QS 
result sets in (for v = 0.01), around A = 1.5. Note that Avar(n) approaches 1/u as 
A — > oo, as predicted by Eq. (34); the absolute value of the difference does not approach 
zero for large A, as it does for the mean. (The relative difference Avar(n)/var(n) does of 
course approach zero in the limit A — > oo.) In Fig. 10 we compare the difference from 
mean-field theory for the third central moment ((n — (n)) 3 ) (mean-field theory yields n 
for this quantity); Fig. 11 is a similar comparison for the fourth central moment [in this 
case the mean-field prediction is (3n + l)n]. The series agrees well with the QS result 
(again, for the case v = 0.01), for A > 2.5. 
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Fig. 11. As in Fig. 9, but for the fourth central moment. 

On the basis of these comparisons we conclude that the perturbation theory developed 
in the previous section is incapable of describing the vicinity of the transition (A ~ 1), but 
agrees extremely well with quasi-stationary properties (including higher moments) above 
a certain, not very large value of A. In this regime, the perturbation approach yields 
accurate analytic expressions for QS properties. 
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1.0 1.2 1.4 1.6 1.8 2.0 

X 

Fig. 12. Correction An to the mean population size versus birth rate A in the MVP with 
v = 0.01. Bold curve: QS distribution; dotted line: lowest order correction in f2-expansion; 
solid line: e series truncated at first order. 

In Fig. 12 we compare the lowest order f2-expansion prediction for An and the corre- 
sponding result of the e series (truncated at order e) with the QS result, for v = 0.01. The 
^-expansion result is seen to be slightly better, although it diverges as A — > 1. Fig. 13 
is a similar comparison for Avar(n). The f2-expansion is again slightly superior. Finally, 
Fig. 14 shows that the fifth-order e series is superior to the second order f2-expansion 
result, as is to be expected. 
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Fig. 13. Correction Avar(n) versus birth rate A in the MVP with v = 0.01. Symbols as 
in Fig. 12. 
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Fig. 14. Correction An to the mean population size versus birth rate A in the MVP with 
v = 0.01. Symbols as in Fig. 12. 



V. MVP WITH INPUT 



In this section we examine the effect of a steady input of individuals at rate 7. For 
7 > the n = state is no longer absorbing, and the system approaches a true stationary 
state as t — > 00. The forward transition rate is now w n>n+ i = 7 + An; the reverse rate 
w n -i,n = n + vn{n — 1) as before. The evolution operator is that of the MV process with 
the addition of a term 7(77 — 1), in the notation of [4]. This corresponds to a new term in 
the action (i.e., the argument of the exponential in Eq. (4)), 

7 f (#' - l)dt' = 7 / ^dt' (49) 
Jo Jo 
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We now proceed as in Sec. II, introducing the shift of variable = n + <p(t). Equating 
the coefficient of ip in the action to zero yields the stationary solution of the macroscopic 
equation, that is, 

n = —(\-l + w) (50) 



where w = a/(A — l) 2 + A*yv. After the shift, the action has the same form as for the MV 
process without a source, but the factors associated with the source and the bifurcation 
are changed. The factor for the source is now: 

—n(X + 1 — w) = uon (51) 

while that for the bifurcation is 1 — w. The factors associated with the conjunction and 
the 4-vertex are —v, as before. 

Consider the contribution to (n^/n due to the lowest order diagram (i.e., the second 
diagram of Fig. 2); this is readily seen to be —uv/w 2 . As before, we may dress the 
diagram (see Fig. 3), which again has the effect of multiplying the above result by the 
factor k — [1 + v jw 2 \~ x . Letting e = 1 — k, the lowest order contribution, when dressed, 
is —uoe. For diagrams with n > 4 both rightmost source and the leftmost conjunction can 
be dressed, leading to the algebraic factor 

. uj s f \ — 1 + & ] 



where m = n — c and u = 1 — w as before. In terms of the parameter e we have 

_rrt 

F alg = (-iy + fhuW m - n {i _ e)m _ 2 (53) 



where 

h = 

X + w-1 



and 



A 2 



u 



V = 



(55) 



With these results the enumerations derived for the original process may be used to 
develop an e series for the MVP with input. The input rate 7 enters via the expression 
for w. Up to third order we find, 
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(^oo) = n 1 — uje + vh \ 2uw — 5— J e 2 
L V wJ 



+ vh -4uw - 10m + 19v + 55m- 



w 



-60(1) 



(56) 



Series is compared with the exact value of (rioo) (from the stationary solution of the 
master equation, obtained numerically), and with the mean- field prediction n in Fig. 15, 
for parameters 7 = « = 0.01. Excellent agreement between series and the exact result 
is found for A > 1.55 and again (see inset) for A < 0.75. In the regions of agreement 
e is small (e < 0.14 for A < 0.75; e < 0.032 for A > 1.55), but it becomes large in the 
intermediate region, attaining a maximum of about 0.96 for A = 1. In this region the 
series prediction deviates strongly from the true value, and can become negative. Similar 
results are found for other values of 7 and v. 




Fig. 15. Stationary mean population size versus birth rate A in the MVP with v = 0.01 
and input at rate 7 = 0.01. Symbols as in Fig. 4. 



VI. CONCLUSIONS 



We apply the path-integral based perturbation theory for Markovain birth-and-death 
processes [3,4] to the Malthus-Verhulst process and a variant of this process that includes 
particle input. At zeroth order the formalism yields a Poisson distribution whose mean is 
given by the mean-field or macroscopic equation. Computational enumeration of diagrams 
allows us to derive the series coefficients for the first four stationary moments of the 
process. The expansion parameter, e, is small when the mean population size (n) is large, 
but is of order unity when (n) is small. In the latter regime the expansion fails, but in 
the former its predictions are in near-perfect agreement with numerical results, despite 
evidence that the series is divergent. Truncating the series at third order, we find excellent 
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agreement with numerical results when the mean population (n) > 35. Our analysis yields 
asymptotic expressions (valid for large population size) for the moments of MVP. 

Our study of the MVP with input shows that the poor behavior of the series, when e 
is not small, is not due the presence of an absorbing state, since input eliminates such a 
state from the process. Since the present approach is based on systematic approximations 
to mean-field theory, one should not expect it to yield useful results where the latter 
is seriously in error, as is the case for A < 1 in the MVP with or without input. For 
parameter values such that the mean-field solution is reasonable, the series provides useful 
corrections to it. Development of a globally accurate perturbation method remains as an 
open challenge, one we hope to explore in future work. 

Comparison with van Kampen's f2-expansion shows that (for the problems considered 
here), the latter method, to lowest order, yields results that are very similar, though 
slightly superior to our first-order expansion. Including higher order terms, the quality 
of the e series improves. In the present case (and in other problems likely to arise in 
stochastic modeling), deriving higher-order corrections is simpler using the diagrammatic 
approach, making our method an attractive alternative to the f2-expansion. 
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